#include "stdafx.h"

#include "PCASolver.h"

#include <iostream>

using namespace Eigen;

PCASolver::PCASolver(Eigen::MatrixXd * data, Eigen::MatrixXd * PC) {
    if (PC != NULL) {
        this->PC = new MatrixXd(*PC);
        return;
    }

    data = new Eigen::MatrixXd(*data);

    int m = data->rows();
    int n = data->cols();

    // substract mean from each dimension
    for (int i = 0; i < m; ++i) {
        double mean = data->row(i).mean();
        data->row(i).array() -= mean;
    }

    // calculate covariance matrix
    MatrixXd covariance = data->matrix() * data->transpose() * ((double) 1 / (n - 1));

    // calculate eigenvectors and eigenvalues
    SelfAdjointEigenSolver<MatrixXd> eigensolver(covariance);
    if (eigensolver.info() != Success) abort();

    VectorXd v = eigensolver.eigenvalues();
    MatrixXd pc = eigensolver.eigenvectors();

    // sort variances in decresing order
    for (int i = v.rows(), j = 0; i > 1; --i, ++j) {
        int maxIndex;
        double maxCoeff = v.tail(i).maxCoeff(&maxIndex);
        v(j + maxIndex) = v(j);
        v(j) = maxCoeff;
        pc.col(j).swap(pc.col(j + maxIndex));
    }

    this->PC = new MatrixXd(pc);
    delete data;
}

PCASolver::~PCASolver() {
    if (PC != NULL) {
        delete PC;
    }
}

MatrixXd * PCASolver::performDimensionReduction(Eigen::MatrixXd * data, int dimensions) {
    return new MatrixXd((PC->transpose() * *data).leftCols(dimensions));
}
